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Abstract 

An effective procedure to determine the optimal parameters appearing in artificial 
flockings is proposed in terms of optimization problems. We numerically examine 
genetic algorithms (GAs) to determine the optimal set of parameters such as the 
weights for three essential interactions in BOIDS by Reynolds (1987) under 'zero- 
collision' and 'no-breaking- up' constraints. As a fitness function (the energy func- 
tion) to be maximized by the GA, we choose the so-called the 7-value of anisotropy 
which can be observed empirically in typical flocks of starling. We confirm that the 
GA successfully finds the solution having a large 7-value leading-up to a strong 
anisotropy. The numerical experience shows that the procedure might enable us to 
make more realistic and efficient artificial flocking of starling even in our personal 
computers. We also evaluate two distinct types of interactions in agents, namely, 
metric and topological definitions of interactions. We confirmed that the topological 
definition can explain the empirical evidence much better than the metric definition 
does. 
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1 Introduction 

Collective behaviour of interacting agents such as flying birds, moving insects 
or swimming fishes shows highly non-trivial properties. We sometimes find 
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a kind of 'beauty' in the quite counter- intuitive and fascinating phenomena 
[1,2,3]. If one wishes to deal with these ingredients by mathematically rigorous 
approach, we sometimes regard each of them as a simple 'particle' without 
size and any specific shape. As a typical example of such 'massive' interacting 
particle systems, a critical phenomenon of order-disorder phase transitions 
with 'spontaneous symmetry breaking' in spatial structures of the so-called 
ferromagnetic Ising system has attracted much attention of physicists. Up to 
now, a huge number of numerical and analytical studies in order to figure it 
out have been done by theoretical physicists and mathematicians [4] . 

On the other hand, for the mathematical modelling of many-particle systems 
having interacting intelligent agents (animals), we also use some probabilistic 
models. For instance, in the research field of physics, Vicsek et.al. [5] proposed 
a flocking dynamics having a simple rule, namely, a given particle (agent) 
driven with a constant absolute velocity at each time step assumes the average 
direction of motion of the particles (agents) in its neighbourhood of radius r 
with some random perturbation added. 

In engineering, a simplest and effective algorithm called BOIDS [6,7] has been 
widely used not only in the field of computer graphics but also in various 
other research fields including ethology, control theory and so on. The BOIDS 
simulates the collective behaviour of animal flocks by taking into account only 
a few simple rules for each interacting intelligent agent. 

Recently, quite a lot of useful flocking algorithms inspired by the BOIDS were 
proposed by a combination of a velocity cooperation with a local potential- 
driven field (for instance, see [8,9]). Among these studies, Olfati-Saber [10] 
provided a remarkable framework for designing of scalable flocking algorithms. 
His framework has three essential factors in the algorithm. The first one is the 
same three essential rules as those in the BOIDS we mentioned just above. 
The second factor is the ability of avoidance of unexpected obstacles appearing 
on the path of flock's movement. The third and the most remarkable one is 
the ability for causing the flock to track the path of a single virtual leader by 
introducing a navigational feedback forth to each agent, namely, all agents in 
the flock are moving according to the information about the virtual leader. 
However, up to now, nobody knows whether such a virtual leader actually 
exists in real flocks or not. 

Hence, it is very hard task for us to evaluate these modelings and also very 
difficult to judge whether it behaves like realistic or not due to a lack of enough 
empirical findings to be compared [11,12,13]. 

As we know from the above issue (doubt) as an example, one of the serious 
problems in studies of any artificial flocking (algorithm) is apparently a lack of 
empirical data to check the validity. Actually, there are few studies to compare 
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the results of the flocking simulations with the empirical data. Therefore, the 
following essential queries still have been left unsolved; 

• What is a criterion to determine to what extent the flocks seem to be real- 
istic? 

• Is there any quantity (statistics) to measure the quality of the artificial 
flocks? 

• Is it possible for us to construct the mathematically defined 'optimal' BOIDS 
in computers? If possible, how does one design the optimal BOIDS in terms 
of some maximization (or minimization) principle of appropriate fitness 
functions? 

From the view point of 'engineering', the above first two queries are somewhat 
not essential because their main goal is to build-up a useful algorithm based 
on the collective behaviour of agents. However, from the natural science view 
points, the difference between empirical evidence and the result of the simu- 
lation is the most important issue and the consistency is a guide to judge the 
validity of the computer modelling and simulation. On the other hand, the 
third query is very important for engineering to solve important problems in 
the real world by using the knowledge of such outstanding abilities of these 
intelligent flockings. 

Recently, Ballerini et. al. [14] succeeded in obtaining the data for such collec- 
tive animal behaviour, namely, empirical data of starling flocks containing up 
to a few thousands members. They also pointed out that the angular density 
of the nearest neighbours in the flocks is not uniform but apparently biased 
(it is weaken) along the direction of the flock's motion. 

With their empirical data by hand, in the previous paper [16], we examined 
the possibility of the BOIDS simulations to reproduce this anisotropy and we 
also investigated numerically the condition on which the anisotropy emerges. 

However, in our previous studies, we checked only the existence of the anisotropy 
and did not check extensively the strength of the anisotropy which is measured 
by the 7-value. To make the matter worse, due to some technical limitations 
of computer simulations, we could not evaluate the n-th order nearest neigh- 
bouring dependence of the 7-value precisely. Namely, due to the following 
three bottlenecks, we could not check the validity of the BOIDS simulations 
by means of the anisotropy measurements: 

• Bottleneck 1 : It is very hard to check whether the result comes from the 
nature of mathematical modelling or from the choice of parameters appearing 
in the model. 

• Bottleneck 2 : It is very difficult for us to generate aggregations with a 
high 7-value which have been observed in a lot of empirical findings. 

• Bottleneck 3 : It is very difficult to evaluate the 7-value precisely for the 
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n-th order nearest neighbour due to the so-called border bias. 

In this paper, in order to overcome the above Bottleneck 1 and Bottleneck 

2, we propose and examine a genetic algorithm (GA) to maximize the 7- value 
which is implicitly regarded as a fitness function of the weights of essential 
three interactions, namely, Cohesion, Alignment and Separation, appearing in 
the BOIDS algorithm. By finding the optimal weights for the BOIDS, we ex- 
pect that the 7-value for the 'optimal' BOIDS is enhanced by the appropriate 
choice of the weights. For the Bottleneck 3, we propose a border-bias free 
procedure to evaluate the 7-value in the computer simulations. 

This paper is organized as follows. In the next section 2, we explain the con- 
cept of anisotropy measurement. The anisotropy distribution map and the 
measurement through the 7-value are also explained. In the next section 3, we 
explain the essential three interactions in the BOIDS. In section 4 and 5, we 
provide the set-up of scale-lengths and time-scale in our computer simulations. 
In section 6, we show the result without optimization as a preliminary. In the 
next section 7, we mention that the selecting the interactions in the BOIDS 
is formulated as an optimization problem to maximize the 7-value as the cost 
function. In section 8, we explain the genetic algorithm to maximize the cost 
function and why we use the algorithm. The results are shown in the next 
section 9. In these modelling, we assume that each agent interacts with the 
other mates within a fixed range of the visual field. In this sense, the model 
should be referred to as metric model On the other hand, one can consider 
the topological model in which each agent interacts with a fixed number of the 
mates. In section 10, we apply our procedure to design the optimal BOIDS for 
the topological model and compare the result with that of the metric model. 
We find that the topological model can reproduce the empirical finding much 
better than the metric model does. In section 11, we discuss the results and 
the last section is summary. 



2 Anisotropy in real and artificial flockings 



In this section, we explain the concept of anisotropy in flockings originally 
proposed by Ballerini et. al. [14] to evaluate the empirical data of starling 
flockings. For the emergence of the anisotropy, we evaluate the strength of the 
anisotropy by the measurement, that is, the 7-value. 
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Fig. 1. Angular distribution map simulated by BOIDS-based modelling (from 
Makiguchi and Inoue (2010) [16]). 

2. 1 Emergence of anisotropy 

Ballerini et. al. [14] measured each bird's position in the flocks of starling 
(Sturnus vulgaris) for 8 seconds in three dimension. To get such three dimen- 
sional data, they used 'Stereo Matching' which reconstructs three dimensional 
object from a set of stereo photographs. From these data, they calculated the 
angle between the direction of nearest neighbours and the direction of the 
flock's motion for all birds in the flock. They measured the angles (0, a), 
where (/) means the latitude (e [—90°, 90°]) of nearest neighbour for each bird 
measured from the direction of the flock's motion, whereas the vertical axis a 
denotes longitude (e [—180°, 180°]) which specifies the position of the nearest 
neighbour for each bird around the flock's motion, of the nearest neighbour 
for all birds in the flock, and plot these angles in the two-dimensional map 
using the so-called Mollweide projection. To put it briefly, this map means the 
density of angular distribution of the nearest neighbour. Inspired by their em- 
pirical findings, we simulate the distribution map by BOIDS simulation [16]. 
The resultant angular distribution map is shown in Fig. 1. This figure clearly 
shows that the density is not uniform but obviously biased. 

2.2 Anisotropy measurement: Formula and generic properties 

To evaluate the strength of the anisotropy, we use the 7-value which is an 
anisotropy measurement introduced by Ballerini et. al. [14]. In following, we 
briefly explain how to compute it and mention the general properties. 

As a matter of convenience, let us first define the vector: 




(i) 
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in the Dime's bracket representation in quantum mechanics [17] as a three- 
dimensional unit vector pointing to the n-th nearest neighbouring agent (bird) 
from an arbitrary argent i. 

We should keep in mind that t appearing in the shoulder of matrix here such 
as A f stands for the 'transpose'. Then, we have the following 3x3 projection 
matrix in terms of the Dirac's bracket: 

M(n) = ^E(l^ (n) )^ (n) l) ( 2 ) 

iV i=i 

whose components are given by 



1 N 

(^ (n) )^=^E(l«l n) ))«(l4 n) )V (3) 

iV i=i 

for a, /3 — x, y, z ) where N stands for the total number of agents in the 
flock. For the above M^ n \ we immediately obtain the normalized eigenvectors 
\U^), k = 1, 2, 3, and one can rewrite the matrix in terms of these bases 

as 



M^ = Y.Mut ] )(ut ] \ (4) 

k=l 

where A fc , k = 1, 2, 3 stand for the eigenvalues of the matrix, that is, M^\U^) = 
Xk\uj^), k = 1, 2, 3, and of course, the rectangular condition (uf^\U®) = S^i 
is satisfied. It should be noted that when the projection matrix becomes 
irregular, we cancel the observation in our calculations of the 7-value. There- 
fore, the probability Pj~ that an arbitrary agent exists in the direction of the 
vector \uj^) is explicitly given by 



P k = \{U^\M^\U^)\ 2 = \l fc = 1,2,3. (5) 

When we put these eigenvalues in a particular order, say, Ai < A 2 < A 3 (this 
reads Pi < P 2 < P3), the vector \U^) is the direction in which the agents are 
more likely to exist, whereas there are fewest agents in the direction of vector 
\u[ n) ). Therefore, if we define the emergence of anisotropy as the absence of the 
birds along the direction of the flock' s motion, the strength of the anisotropy 
is naturally measured by the inner-product of the vector pointing to the flock's 
movement and the eigenvector having the lowest eigenvalue. 

To evaluate the anisotropy measurement more explicitly, let us define the 
eigenvector of the lowest eigenstate by |Vt^ n) ), that is, {W^ n) \M {n) \W^ n) ) = 
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min& Afc. Then, the anisotropic measurement is calculated from the | W^) and 
the vector \ V) that points to the direction of the flock's movement (the velocity 
of the center of mass in the flocking \V) = (1/AT) Yh=i Vi) as 



Obviously, the above j t is dependent on the time t through the time-dependence 
of \ V) and Thus, we define the anisotropic measurement 7 by averaging 

over the 'observation-time' with the infinite length T —> 00 as 



As it is impossible to take the infinite observation-time limit T 00 in 
computer, we replace the limit by finite observation-time, say, T = 80 for the 
upper bound in the sum (7). 

We should also keep in mind that the 7- value depends on the choice of initial 
conditions and one should take the average over the distribution of initial 
conditions. However, we might expect that the 7- value calculated for a single 
realization of initial conditions, say, 7 is identical to its average E^. [• • • ] over 
the initial condition, namely, 7 = E^.py] in the limit of TV — >• 00. As we shall 
show later, the number of agents in our simulation is too small N = 100 <C 00 
to satisfy the above condition. Therefore, we calculate the average of the 7- 
value for 1000-independent initial conditions. 

It should be noted that the 7-value takes any positive values in the range 
< 7 < 1. For the case of 7 = 0, the nearest neighbour is more likely to exist in 
the direction of flock's movement, namely, (W^IV) = 0. On the other hand, 
the 7=1 implies that the nearest neighbour exist in the two-dimensional 
plane which is perpendicular to the flock's movement with probability 1, that 
is to say, \W^ n) ) = ±\V) resulting in \(W^ n) \V)\ 2 = 1. We should notice that 
for the eigenvector \S^) and {S^) having the largest and the second largest 



eigenvalues, \ (S^\V)\ 2 +\(Si n) \V)\ 2 = 1 for 7 = and (S^\V) = (Si n) \V) = 



for 7 = 1 should be satisfied. 

After simple algebra, one can show that the 7 takes 1/3 when there is no 
spatial bias in the direction of the nearest neighbours (namely isotropy). This 
means that the anisotropy emerges when the following condition is satisfied. 



7t = \(W( n) \V)\ 2 . 



(6) 



7 = E t [ 7t ] = lim ~Y.lt- 




(7) 



_ 1 

T ^ Tuniform = 77 



(8) 
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In fact, the above inequality is easily confirmed. The 7- value for the uniform 
distribution of the position (0, a) ) where (j) and a are the same variables of 
angles defined in the previous subsection, for a given vector \V) ) namely, the 
7-value for a) = (47r) _1 is easily calculated as 

Tuniform = / p(0, *) # da \ {W™ \ V ) | 2 

J sphere 

_ 1 

~ 3 

where we used |(Vl^^|y)| 2 = cos 2 (j) cos 2 a. Therefore, the distribution of the 
n-th nearest neighbours has an anisotropic structure when the 7-value is larger 
than 7 un iform = 1/3, namely the condition for the emergence of the anisotropy 
is explicitly written by (8). 

Ballerini et. al [14] measured the 7-value up to the n-th order of nearest 
neighbours for two kinds of empirical data having different numbers of agents 
and show the n-dependence of the 7-value. From their plot, we clearly find 
that 7-value takes larger than 0.8 for n = 1 and the value remains larger than 

Tuniform 

= 1/3 up to n = 6. In our previous studies [16], the n-dependence 
of the 7-value was evaluated for the data generated artificially from BOIDS 
simulations, however, we could not overcome the problem of border bias which 
was mentioned as Bottleneck 3 in the previous section. In this paper, we 
introduce a way to overcome this technical difficulty and attempt to measure 
the 7-value more precisely. 



3 Essential three interactions in BOIDS 



To make flock simulations in computer, we use the so-called BOIDS which 
was originally designed by Reynolds in 1987 [6]. The BOIDS is one of the 
well-known mathematical (probabilistic) models in the research fields of CG 
and animation. Actually, the BOIDS can simulate very complicated animal 
flocks or schools although it consists of just only three simple interactions for 
each agent in the aggregation: 

(c) Cohesion: Making each agent's position Xi (i — 1, • • • , N) toward the 

average position of neighbouring flock mates, 
(a) Alignment: Keeping the velocity of each agent Vi(i = 1, • • • , N) the 

average value of neighbouring flock mates, 
(s) Separation: Making a vector of each agent's position Xi (i — 1, • • • , N) 

to avoid the collision with the neighbouring flock mates. 



1 r 2 s 

= — / cos a da cos 6 dS 

47T J-7T J-l 



(9) 
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Each agent decides her (or his) next direction of migration by compounding 
these three vectors of interaction. In addition to this, it is important for us to 
bear in mind that 'local flock mates' mentioned above denotes the neighbours 
within the range of view for each agent. We explain this view and other settings 
of our simulation in the next section. 



3.1 BOWS dynamics 



For simplicity, we define 'neighbouring mates' by all metes which exist within 
the visual field with a radius R and for each mate categorized as the neigh- 
bouring mates, we calculate the interactions of Cohesion and Alignment as 
normalized unit vectors. In addition, we evaluate the interaction of Separa- 
tion as a unit vector pointing to the direction of fading-out from the mates. 
Each agent i updates its own velocity vector V{ and the position Xi by the 
following recursion relations. 



V t {l + l) = vfe^{l) (10) 
X i (l + l) = X i (l) + Vi(l + 1) (11) 

where I denotes time step in our simulations and we discretized the infinites- 
imal time as a unit time step Al = I + 1 — I = 1 in the definition of velocity 
V % = dX z /dl ~ {Xi(l + Al) - Xi(l)}/Al = Xi(l + 1) - Xi(l) = Vi(l + 1) to 
obtain (11). e^{l) denotes a unit vector pointing to the direction to which the 
agent i should move according to the BOIDS. The e^(l) is explicitly given by 



\j 1 v ( ( ?(i)+j 2 vf(i)+j 3 v ( s i) (i)\ +v \Vi(i)\ 



Jiv$(i)+j2v ( ?(i)+j3v%\i) , Vm 
|Jivg ) (o+J 2 v«(o+J 3 'y^ ) (OI '\Vi(i)\ 



(12) 



with 



4°(0 = 
^(0 = 

4°(0= 



^ =1 9(R-r aj )J 3 (l) 



Xi(l) 



0(fl-ry).X:,(!) 



Ef =1 6(^-^)^.(0 



-Xi(l) 



£? =1 0(^-^)^.(01 



Ef =1 0(R 



(13) 

(14) 
(15) 
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where we defined as the square distance between agent i and j as 



r y = \Xi(!) - Xj(l)\ = ^{X t (l)-X,(l)r. (16) 

0(- • • ) denotes a step function. Therefore, Ylf=i Q(R — Tij) stands for the 
number of 'neighbouring mates' for the agent i and the number is obviously 
dependent on the agent i. 

A balance parameter 77 appearing in (12) determines the weights of two dis- 
tinct modifications e^{l) for the velocity vector for the agent z, namely, the 
'BO IDS- driven' correction term ~ J\V^q J^v^ [1) + J^v^\l) and the vector 
conservation term ~ Vi(l). Hence, the vector Vi{l) is conserved for 77 3> 1, 
whereas the dynamics of Vi{l) becomes purely BOIDS-driven for <C 1. 
Therefore, the choice of r] is regarded as a kind of 'inert effect' in the dynam- 
ics of BOIDS. The value itself should be determined empirically. However, due 
to the lack of such useful information about the inertia in real flockings, here 
we simply set 77 = 2 by ad-hoc manner. 

From the above definition of (13), we easily find that Vq\1) = —v$\l) and one 
of these two distinct effects is completely cancelled in the BOIDS dynamics 
(10)(11) as ~ any choice of J u J 3 . To correct this undesirable 

situation, we slightly modify the v$\l) as follows. 

p (o m e(i?o-r t7M )(x^(0-x t (0) 

SU mRo-r^XX^-Xm 1 ) 

where j(n : i) denotes the n-th nearest neighbouring mate of the agent z and 
it is explicitly given by : 



j(n : z) = axgrnax^^zj^^.. r %3 (18) 
with j(0 : i) = i. Then, j(l : i) means the nearest neighbouring site 



j(l : z) = argminr^ (19) 

3 

for each z. Namely, the separation v$ (I) acts if and only if the distance between 
the agent z and the nearest neighbouring mate j(l : i) is lower than the radius 
of separation range R . 

The vectors Vq\1),v^(1) and v$\l) denote the components caused by the 
interactions Cohesion, Alignment and Separation for agent i at time /, respec- 
tively. We should keep in mind that \v§(l)\ = \v^(l)\ = \v%\l)\ = 1 holds 
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and J = ( Ji, J2, J3) stands for the set of weights for three interactions, namely, 
Cohesion, Alignment and Separation. We should keep in mind that from (13), 
(14) and (15), the weight J 2 is a 'dimension-less' variable, however, J\ and J 3 
have inverse-time dimension ~ (time step) -1 . 

We should keep in mind that the above definition of j(n : z), the 7- value is 
calculated by setting \u^) = u[ n) = (V^ - V^V^ - V x \ in (3)(6) 
and (7). 

From equation (10), we are confirmed that the amplitude of velocity vector of 
agent i at time / + 1 is identical to the average amplitude of velocity vectors 
for neighbouring mates in the previous time step I as 



K0 _ S ga 6(^-^)1^( 01 
Ef =1 e(R - r y ) 



\VS + 1)| = V\ l> ee ^^7^ V ■ (2°) 



The above update rules (10), (11), (12), (13), (14) and (17) are our basic 
dynamical equations to be evaluated numerically. 

Obviously, the behaviour of the artificial flockings strongly depends on the 
choice of the weights, however, there is no extensive study to investigate to 
what extent the behaviour changes quantitatively by changing the weights. 
From the fact in mind, in this paper, we propose a systematic algorithm to 
determine the weight by using the evolutionary computation such as GAs to 
maximize the 7-value as a fitness function. 



4 Scale-lengths in BOIDS simulations 



We here explain how we set several scale-lengths appearing in our simulations. 
In our previous studies, we determined them without any justification from 
the empirical evidence, however, here we attempt choose the scale lengths by 
taking into account the data from the reference [15] in order to realize the 
artificial flockings as realistic as possible (We summarize these variables in 
Table 1). 

However, some parameters is not determined by empirical data, for instance, 
the range of interaction, the Frame-Rate (FR) [fps:frame-per-second] and so 
on. Therefor we set these parameter by the subjective view point and some 
regards for the calculation cost, for example, the number of agents TV = 100, 
the radius of the visual field R = 3 x R where R = 1.09 [m] denotes the 
radius of separation range and the FR, which will be explained in the next 
section in detail, is 200 [Hz]. 
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A set of scale-lengths in our BOIDS 



Number of agents (N) 


100 


Body-Length (BL) 


0.2 [m] 


Wmg-bpan (Wo) 


r\ a r 1 

(J. 4 [mj 


Radius of Separation Range (Rq) 


1.09 [m] 


Radius of Visual Field (R) 


3 x Rq [m] 


Initial Speed Average {V') 


10.10 [m/s] 


Initial Density of the Aggregation (p) 


0.13 [m" 3 ] 



Table 1 

A set of scale-lengths in our flock simulation. Variables other than the Number of 
agents and the Radius of Visual Field are based on empirical data by Ballerini et. 
al (Event 29-03 in Table 1 of [15]). 

5 On the time-scale in BOIDS simulations 



In our previous study [16], we defined the unit time (frame) by 0.1 [sec]. In 
this paper, we shall define the frame based on the so-called Frame-Rate(Fi? = 
200 [Hz]). In order to consider the consistency with the empirical data analysis 
by Ballerini et. al. [14] in which they used 0.1 [sec] for a unit frame, we evaluate 
the 7- value every FR/ 20 frames and the distance covered by each agent per 
frame, that is, the average of flock's velocity V is also determined from the 
empirical evidence of velocity V' [m/sec] as V = V' /FR [frame -1 ]. 

In our previous work [16], we also chose the initial velocity for each agent 
from a uniform distribution having a finite support. However, this procedure 
might cause some difficulties, namely, we might encounter the 'breaking-up' 
of the flocking to several small groups due to synchronization in their speeds 
of convergence. In general, it is very difficult for us to control the speed of the 
flocking (the speed of the center of mass) after each agent's speed converges 
when we determine the initial speed of each agent by a random number from 
a uniform distribution. To overcome this type of difficulties, we sample the 
initial value of each agent's velocity Vi from the following Gaussian with mean 
V = 10.10 [m/sec] and variance a 2 = (V — l)/3, namely, 



By this setting, we are confirmed that the speed of flocking actually converges 

toy'. 



(Vi - v f 

2a 2 



(21) 
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Crowded 


(1,0,0) 


0.332 [0.0816] 


100% 


Synchronized 


(0,1,0) 


0.319 [0.292] 


17.93% 


Spread 


(0,0,1) 


0.347 [0.304] 


0% 



TablG 2 

Resulting 7- value and the SD (Standard Deviation), and the FC (Frequency of 
Collisions) are shown for three different ad-hoc choices of the weights (Ji, J2, J3). 

6 A preliminary: Simulations without GA 

In this section, we show the results without any searching of the optimal 
weights of interactions by genetic algorithms as a preliminary. 

6.1 Preliminary results 

In the above setting of the problem, we attempt to evaluate the 7-value using 
the same way as our previous study [16]. We control the weights of three 
interactions in the BOIDS J = (Ji, J2, J3) to generate typical three cases, 
namely, Crowded for (^1,^2,^3) = (1,0,0), Synchronized for (^1,^2,^3) = 
(0,1,0) and Spread for (Ji, J2, J3) = (0,0,1). We list the 7-values and the 
corresponding frequency of collisions (FC) in Table 2. We also show the angular 
distribution maps in Fig. 2. From these table and figure, we clearly find that 
in all cases, the anisotropy is not observed at all. 

We next choose the weights J\ and J 2 as J\ = 1, J 2 = 5 which we used in the 
previous study [16] as an appropriate choice to produce the anisotropy, and we 
shall vary the J 3 from 0.2 to 2.0 to evaluate the 7-value and the corresponding 
frequency of collisions as a function of J 3 . We plot them in Fig. 3. From this 
figure, we find that the frequency of collisions decreases monotonically as J 3 
increases and it converges to zero around J 3 = 1.4. On the other hand, the 
7-value takes its maximum 0.70 around J 3 = 1.2. From these observations, we 
conclude that we should choose the weight as J = ( Ji, J 2 , ^3) = (1, 5, 1.4) and 
whose 7-value will be about 0.70 . 



7 Optimization to design artificial flockings 

In the previous section, we show on what condition the anisotropy emerges 
by setting the weights for three essential interactions by ad-hoc manner as a 
preliminary. Here we mention that the procedure to determine the interactions 
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Fig. 2. From the top to the bottom, the angular distribution maps for Crowded, 
Synchronized and Spread cases. In all cases, the anisotropy is not observed at all. 

can be regarded as a kind of optimization problems. 



7. 1 Optimization under two essential constraints 



In real flockings, it might be a very serious problem for each agent how to 
refuse collisions with the other mates during the flock is moving. In addition, 
it is very hard for us to say that the flock splitting into several sub-flocks is 
also 'realistic' flocking. Therefore, we should design the BOIDS simulations so 
as to avoid these two unexpected accidents, namely, 'collision' and 'breaking- 
up'. To realize the simulation in which there are no collision and breaking-up, 
we introduce two constrains into the optimization problem. 



We first define the 'collision' as the case in which the distance between an 
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Fig. 3. The 7-value and the corresponding frequency of collisions as a function of J3 
we set J\ = 1, J2 = 5. The frequency of collisions drops to zero around J3 = 1.4. This 
fact makes us determine to choose the appropriate set as (Ji, J 2 , J3) = (1, 5, 1.4). 



arbitrary agent z and its 1-st nearest neighbouring mate j(l : z), say, ^j^y is 
shorter than their body length BL= 0.2. We also assume that the 'breaking- 
up' occurs when the distance between the center of mass and the most far 
agent from it becomes longer than a given constant. It is naturally imagined 
that taking into account the 'collision' and making the algorithm to avoid it 
are essential issues not only for artificial flockings but also for flockings in the 
real world. It should be noted that the conventional flocking simulations based 
on the 'particle models' (see for example [10]) in which the size of the mate is 
neglected cannot deal with the 'collision'. We are also confirmed that avoiding 
the 'breaking-up' also might be an essential factor to decide the size of the 
flocking. 

Thus, we might use the 7-value as a cost function (energy function) to be 
minimized under 'zero-collision' and 'zero-breaking-up' constraint to deter- 
mine the three essential interactions J = ( Ji, J2, Js)- Namely, we should solve 
the following optimization problem with the cost. 



^(J) = 7 (J) + A 1 A^(J) + A 2 i3(J), Ai,A 2 ->oo (22) 

where we defined N and B as the number of the collisions and breaking-up, 
respectively. The Ai, A 2 stand for the Lagrange multipliers. In other words, the 
optimal interactions J opt is given by 
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J opt = argmax j lim E(J). 



(23) 



Since Reynolds proposed the BOIDS, quite a lot of the modifications or the 
variants were constructed in terms of engineering, however, no studies con- 
cerning the systematic determination of the essential three interactions in the 
algorithm from the view point of empirical observation on real flockings such 
as starlings. Therefore, here we formulate the procedure to determine the in- 
teractions as optimization problems having the 7- value as the cost under two 
essential constraints. To solve the optimization problem by means of the con- 
ventional tools, say, genetic algorithm (GA for short), we might design the 
BOIDS more systematically. In the next section, we shall examine the GA to 
solve our optimization problem to determine the optimal set of the interactions 
in the BOIDS. 



8 Genetic algorithms 

Here we apply the GA to the determination of the weights of the interactions 
J in BOIDS simulations. Before we show the results, we shall briefly explain 
the motivation to use the GA and the outline of the set-up and the procedure. 
The details of the GA shall be explained in Appendix A. 

8. 1 Why do we use the GA ? 

The GA is a stochastic method to obtain a candidate of the solution having 
the highest possible fitness in the complicated fitness function with multi- 
valley structures. In GAs, one codes the candidates of the solution by a set of 
vectors, each of which is referred to as a 'gene configuration' (a genetic code). 
Then, we make several operations, namely, Crossover, Mutation and Selection 
to regenerate gene configurations having relatively high fitness values [18]. 

As a study to determine the weights for the interactions in the BOIDS, Chen 
et. al [19] proposed the so-called Interactive genetic algorithm (IGA). How- 
ever, we should mention here that they used the fitness function which is 
constructed subjectively, and in this sense, their approach is essentially differ- 
ent from ours. This is because as we already mentioned, we use the 7- value 
which is a measurement introduced by empirical findings [14]. 

Of course, there are a lot of optimization methods and we do not have to use 
the GA to obtain the solution to maximize the 7-value. However, we might 
assume that the agent (bird for instance) acquired such an intelligent way to 
behave as 'flock' during their process of evolution and this assumption makes 
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us use the GA. The justification of using the GA is very difficult to show 
and it might be impossible to prove the validity of the above assumption 
theoretically. Nevertheless, here we use the GA as a first attempt to design 
the optimal BOIDS based on the maximization principle of the 7-value as a 
fitness function. 



8.2 Procedure of the GA 



In this subsection, we shall explain the outline of the procedure of the GA. 
In our GA, we use the three weights of the BOIDS, namely, J = (Ji, J 2 , ^3) 
as gene configurations. Each of the components Ji denotes a 'chromosome' or 
simply 'gene' and takes the value in the range [0.001, 0.999] and the minimum 
value of changing the state is set to 0.001, namely, we here vary the value of 
each component by Ji — >• ±0.001, i — 1, 2, 3 for each step of three operations 
we mentioned above. 

To operate the Selection, we need the value of the fitness function, namely, 
the 7-value for the nearest neighbouring agent (7-value defined by (7) with 
n = 1). To evaluate it, we use the time-averaged 7-value K t [y t ] (see (7) for its 
definition) which is calculated by sampling positions of the mates every 0.1 
[sec] during 8 [sec] (T = 80 data points are needed to evaluate the 7-value for 
each update of the gene configurations). The details of the total procedure of 
the GA is given in Appendix A. 

We also explain a border-bias free (the 'border-bias' was already mentioned 
in the section of introduction as Bottleneck 3) procedure to evaluate the 
7-value in the computer simulations in Appendix B. 



9 Results 



In Table. 3, We show the highest 7- values for three independent runs, which 
are referred to as Case 1,2 and Case 5, and corresponding weights for the 
interactions. It should be noted that we normalized the weights J = ( Ji, J 2) J3) 
so as to make the maximum Ji among the three i = 1,2,3 unity. For each 
run, the highest 7-value is larger than 0.79 and corresponding weights J = 
(Ji, J2^Js) take similar values for all cases. In following, we investigate the 
optimization process for Case 3. 
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Uase 


7- value 


7 7 7 


Case 1 
Case 2 
Case 3 


0.795 
0.797 
0.797 


0.270 0.640 1 
0.234 0.699 1 
0.190 0.895 1 



Table 3 

Resultant sets of weights for three interactions. 
1 
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Fig. 4. Evolution of 7-value in generations. 
9. 1 Optimization process of GA 

We show the minimum, average and maximum of the 7-value for each gen- 
eration in Fig. 4. From this figure, we find that these all values converge to 
7 ~ 0.8 from the initial state (7 ~ 0.4). We next show the time evolution of 
the weights J = (Ji, J2, ^3) for the best possible gene configuration having the 
highest 7-value in Fig. 5. From this figure, we find that each weight changes 
its state during the GA dynamics and this result tells us that optimal gene 
configuration can be successfully generated by our GA procedure. We also plot 
the distribution of 7-value at the initial generation (Fig. 6 (left)) and at the 
final generation (Fig. 6 (right)). From these two panels, we are confirmed that 
the gene configurations having relatively high fitness values are generated and 
they actually survive until the final generation. However, the final distribution 
has a finite deviation instead of a single delta peak. This means that we could 
not find the optimal gene configuration with probability 1. 
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Fig. 6. Left: Histogram of 7-value for the initial. Right: Histogram of 7-value for 
the final. We are confirmed that the GA successfully finds the value which is close 
to the highest possible 7-value with a high probability. 
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Fig. 7. Histogram of weight of each interaction for the final gene set. Left: weight 
of 'Cohesion'. Center: weight of 'Alignment'. Right: weight of 'Separation'. 
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Fig. 8. The 7-value as a function of n calculated by the 'optimal' BOIDS simulation. 
The inset stands for the corresponding angular distribution. 

9.2 The 7 -value for the optimal BOIDS 

Here we examine the n-th nearest neighbouring agent's 7-value for the optimal 
BOIDS having the optimal weights obtained in the previous section. We carry 
out 1000-trials to evaluate the 7-value for each n from n = 1 up to n = 25. 
We also show the angular distribution map for n = 1 and check the behaviour 
by graphical user interface (GUI). The results are shown in Fig. 8. From this 
figure, we find that the 7-value for n — 1 takes the highest value which was 
also observed in the empirical data analysis [14]. We also find from the GUI 
that a realistic flocking's behaviour in which the distance between nearest 
neighbouring agents is not zero ('zero-collision') but finite is achieved for the 
BOIDS with optimal weights of the interactions . 



9.3 Anti-anisotropy effect and its possible explanation 

In Fig. 8, we find that the 7-value becomes lower than the isotropic limit 7 = 
1/3 for 7 < n < 14. This 'anti-anisotropy' effect can be explained from the 
view point of geometric structure of the flocking as follows. When we assume 
2-Dimensional Field and an arbitrary agent is surrounded by the other six 
mates as shown in Fig. 9 (this hexagon-shape is made by equilateral triangle 
with neighbours), the fifth and the sixth nearest neighbouring mates are more 
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Fig. 9. A reasonable explanation of the 'anti-anisotropy'. The fifth and the sixth 
nearest neighbouring mates are more likely to exist in the direction of flock's motion 
(left). As the result, the 7- values for the 5, 6-th orders of the nearest neighbouring 
become lower than the isotropic limit 1/3 (right). 

likely to exist in the direction of flock's motion (they are indicated by '5' 
and '6' in Fig. 9 (left)). As the result, the 7-values for the 5,6-th orders of 
the nearest neighbouring become lower than the isotropic limit 1 /3 (see Fig. 9 
(right)). We should keep in mind that the 'anti-anisotropy' effect might appear 
much more clearly for the flocking being longer (in the moving direction) than 
is wide. 

To confirm this assumption much more explicitly, we evaluate the third-power 
of average distance R between an arbitrary agent and the n-th nearest neigh- 
bouring mate as a function of n. We show the result in Fig. 10. From this 
figure, we clearly find that the R 3 is almost constant up to the 8 ^ 10-th near- 
est neighbour. This numerical result tells us that an arbitrary agent might be 
surrounded by 8 ~ 10 mates leading to the regular-polygon structure. 

In the empirical data analysis of starling flocking, such 'anti-anisotropy' has 
never observed. Therefore, we might conclude that it is very hard for us to 
accept the regular-polygon structure around an arbitrary agent although our 
algorithm presented in this paper suggested the possibility (Appendix A). 
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Fig. 10. Third-power of average distance R between an arbitrary agent and the n-th 
nearest neighbouring mate as a function of order of neighbour n. We are confirmed 
that the R 3 is almost constant up to the 8 ~ 10-th nearest neighbour. 

10 Topological definition of neighbours in BOIDS 

In the previous sections, we attempted to construct the BOIDS algorithm in 
which each agent interacts with each other when the distance between them 
is shorter than the constant radius of the visual field R. In this sense, we 
utilized the metric definition of neighbours in the BOIDS. In this definition 
of neighbours, the number of agents who interact with an arbitrary agent is 
not constant but apparently fluctuates. As we mentioned, the resulting 7- value 
shows 'anti-anisotropy' due to the regular-polygon structure around the agent. 
Unfortunately, in real flockings, we have never observed such 'anti-anisotropy' 
so far. This empirical fact tells us that it is less likely to exist such regular- 
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polygon structure in the real flocks. 



In fact, Ballerini et al. [14] suggested that a bird in the real starling flock 
interacts with a fixed number of neighbours (about six or seven neighbours). 
From this empirical findings, we conclude that the neighbours in the flocking 
should not be defined by the metric sense but it should be determined by 
the topological sense. Obviously, the topological definition of the neighbours 
is completely different from the metric definition which was adopted in our 
modelling of artificial flockings. 

Hence, this empirical fact also gives us motivations to reconsider the metric 
definition of the neighbours in the BOIDS, namely, here we assume that the 
wrong definition of the neighbours causes the ' count er-empirical' results in our 
computer simulations. 

In this section, we shall reconstruct our BOIDS algorithm by taking into ac- 
count the above empirical fact, namely, topological definition of the neigh- 
bours. 



10.1 Topological model 

To avoid confusion, we first remind readers of two distinct definitions of neigh- 
bours. 

In Fig. 11, we show the cartoons for these two definitions. The left panel shows 
the metric definition of neighbours which we used in the previous sections. As 
we explained, each agent interacts with the others when the distance between 
mates becomes shorter than the constant radius of the visual field R. In the 
case shown in this panel, the agent located at the center of the circle interacts 
with four neighbours. On the other hand, the same agent as in the left panel 
interacts with six neighbours in the case of the right panel. The definition of the 
neighbours shown in this right panel is referred to as topological. Apparently, 
in the topological definition of neighbours, the number of mates interacting 
with a given arbitrary agent is a fixed constant and we define the number as 
n c . Thus, the n c in the right panel of Fig. 11 is n c — 6. 

From now on, the model constructed by means of the metric definition of 
neighbours is referred to as metric model, whereas we call the model based on 
the topological definition as topological model. 

In our topological modelling, we set the number of interacting mates n c = 6, 
which is suggested by empirical data analysis by Ballerini et. al. [14]. 

In order to construct the effective BOIDS simulation based on the topologi- 
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Fig. 11. Two types of the definition for interacting neighbours. The left panel shows 
the metric definition, whereas the right panel corresponds to the topological defini- 
tion. The number of mates interacting with a given arbitrary agent is n c = 6. 

cal definition of neighbours, we should introduce the following new types of 
interactions into our previous BOIDS. 

(tc) Topological Cohesion: Making a vector of each agent's position Xi (i = 
1, • • • , N) toward the average position of neighbours in the topological 
sense, namely, the average position over the neighbouring mates up to 
the n c -th nearest neighbour. 

(gc) Global Cohesion: Making a vector of each agent's position Xi (i = 
1, • • • , N) toward the center of mass in the flocking in order to prevent 
the flock from splitting into more than two distinct clusters. 

(ta) Topological Alignment: Keeping the velocity of each agent V{ (i = 
1, • • • , N) the average value of topological neighbours (up to the n c -th 
nearest neighbour). 

We also slightly improve Separation as 

(ms) Modified Separation: Making a vector of each agent to avoid the col- 
lision with mates up to the n c -th nearest neighbour in the topological 
sense. 

It is expected that the above Modified Separation enables us to avoid mak- 
ing regular-polygon structures in artificial flockings. Namely, in our previous 
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BOIDS simulations, the regular-polygon structure might be induced by met- 
rically defined Separation which acts for the only 1-st nearest neighbour 
mate. 



10.2 BOIDS dynamics 



Hence, our topological model is described by the following update rules 



V t (l + l) = vfe^(l) 



B 



X i (l + l) = X i (l) + V i (l + l) 



(24) 
(25) 



where e B < (I) denotes a unit vector pointing to the direction to which the agent 
i should move according to the above interactions of BOIDS and explicitly 
given by 
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(28) 
(29) 

(30) 



where i?W denotes the square distance between the agent i and the n c -th near- 
est neighbouring mate. Therefore, the number n c should be defined explicitly 

by 
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N 



"c = £e(/2g-^.) (3i) 

for all i because 0(- • • ) survives only for the j satisfying < and the 
number of such j is just n c from the definition. The balance parameter r] is 
set to the same value 2 as in the case of the metric model. 

From the above definition of (27), we easily find that v^ c (l) = — ^ms(0 an d 
one of these two distinct effects is completely cancelled in the BOIDS dy- 
namics (24) (25) as ~ any choice of J u J 3 . To modify this 
undesirable situation, we slightly change the v$\l) as follows. 



?MsU lE^e^-r.^cx^o-x.co)! 



where j(i : n) is given by the definition (18), and R n means the Separation 
Range for the n-th nearest neighbour mate. 

From the empirical evidence [14], we set R® as 



i?« = |r^|n 1/3 (33) 

where we define R$ = R = 0.73 and is selected as a Gaussian variable 
with mean R and unit variance. 

From equation (24), we should notice that the amplitude of velocity vector of 
agent i at time / + 1 is identical to the average amplitude of velocity vectors 
for the topologically defined neighbouring mates, that is, the mates up to the 
n c -th nearest neighbours in the previous time step I as 



\Vi(l + 1)| =Vf = ^ Ll9 ^ r ^'.n))\ V J^ l )\ ^ ^ 

The above update rules (24), (25), (26), (27), (28), (32) and (30) are our basic 
dynamical equations to be evaluated numerically. 

We list a set of scale-lengths appearing in our simulations in Table 4. In both 
Table 1 (metric model) and Table 4 (topological model), we chose these scale- 
lengths from the empirical data [14]. However, we should keep in mind that 
the choices of R are different in both cases. In the metric model, we used 
Ro = 1.09 which is chosen from Event 29-03 in the reference [15], whereas 
Ro = 0.73 in the topological model comes from Event 28-10 in [15]. 
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A set of scale-lengths in the topological model 



Number of agents (N) 


100 


Body-Length \BL) 


0.2 |mj 


Wing-Span (WS) 


0.4 [m] 


Radius of Separation Range (Rq) 


0.73 [m] 


Initial Speed Average V' 


11.10 [m/s] 


Initial Density of the Aggregation (p) 


0.54 [m" 3 ] 



Table 4 

A set of scale-lengths in the topological model. Variables other than the Number of 
agents are based on empirical data by Ballerini et. al. (Event 28-10 in Table 1 of 
[15]). 



For the topological model obtained by the above modifications, we utilize 
the GA to find the weights of the four interactions ((tc),(gc),(ta) and (ms)). 
Then, we numerically evaluate the 7- value and the the third-power of average 
distance (R 3 ) between an arbitrary agent and the n-th nearest neighbour as 
a function of order of neighbour n. 



10.3 Results 



We show the result in Fig. 12. We plot the 7- value as a function of order n 
(left panel) and the R 3 as a function of order of neighbour n (right panel) for 
n c = 6. 

In the left panel, we find that the anisotropy emerges (7 > 1/3) up to n c = 6 
and 'anti- anisotropy' disappears as we expected. From the right panel, we are 
also confirmed that the regular-polygon structures never emerges in the flock 
because the R 3 monotonically increases as the n increases. Of course, here 
we used the empirical finding R n = R^n 1 / 3 [14] to determine the R$ in the 
v ms(0 f° r our topological modelling, however, the result might justify that 
the flock in the metric model behaves as a 'crystal form', whereas the flock in 
the topological model looks like 'gas' which is much closer to real flockings. 

We also calculate the so-called integrated conditional density T(r) and pair dis- 
tribution function g{r) introduced by Cavagna et. al. [24]. These two quantities 
are explicitly defined as 



27 



1 




T-BOIDS N=1Q0 



Order of the neighbor(n) 



Order of the neighbor(n) 



Fig. 12. The resulting measurements for our topological model. The left panel shows 
the 7- value as a function of order n, whereas the right panel is plotted as third-power 
of average distance R between an arbitrary agent and the n-th nearest neighbouring 
mate as a function of order of neighbour n. We set n c — 6, which is indicated by 
empirical evidence [14]. 



1 N-(r) 



g(r) 



1 1 " c 



Aitr 2 n, 



(35) 
(36) 



where A^(r) denotes the number of points in the sphere with the radius r 
centered in z, and n c stands for the number of individuals in the sphere, 
is the absolute distance between the center i and the neighbour j. Here we 
should notice that these two quantities are related each other and satisfy 



0(r) = r(r) + --^. (37) 

Hence, we directly evaluated T(r) from our simulations and calculated the 
g(r) by means of the above equation. 

We show the results in Fig. 13. From this figure, we find that both T and g for 
the metric model suddenly increase around the radius of view field R = 1.09 
and the behaviour is completely different from the empirical evidence [24] . On 



28 



the other hand, these quantities for the topological model gradually increases 
around r = 0.4<i? = 0.73 and the behaviour is very close to the empirical 
evidence (see Fig. 1 and Fig. 5 in [24]). 
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Fig. 13. The integrated conditional density T(r) (upper panels) and the pair dis- 
tribution function g(r) (lower panels) for the metric model (left panels) and the 
topological model (right panels). 



From these numerical results, we conclude that our topological model recon- 
structs the essential features confirmed from empirical findings of starling. Al- 
though some tiny gaps have been left, for instance, the 7-value of the nearest 
neighbour is lower than the empirical evidence, one can say that our topolog- 
ical model is more likely to be 'realistic' than the metric model. 
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11 Discussion 

In this section, we discuss several issues on our results. 
11 A On the optimal gene configuration 

From Table 3, we find that the order of the strength of the weight Ji, J 2 and 
J 3 is J3 > J 2 > J\ in any run of the GA. The condition J 3 > J\ is needed 
to require the 'zero-collision' constraint, whereas the condition J 3 > J 2 > J\ 
makes the 7- value larger than that for the condition J3 > J±. Actually, in 
our previous study [16], we carried out the simulations under the condition 
J3 > J 1 > ^2 ? which was determined by hand, and found that the 7- value 
for n = 1 is around 0.7 which is apparently lower than the result of empirical 
findings 7 ~ 0.8 [14]. This fact is a reasonable advantage of our approach 
based on the GA to maximize the 7-value. 

On the other hand, for the topological model, we find the order of interactions 
as J2 > J3 > J4 > J\ as shown in Fig. 14. 

11.2 Theoretical upper bound of the 7 -value 

In our simulations, we used the 7 value as a 'cost' for the optimization prob- 
lem to determine the three interactions in BOIDS. This procedure is possible 
because the 7-value itself is described as a function of these interactions im- 
plicitly, namely, 7 = 7(Ji, J2, ^3)- Thus, by using the GA (the choice of GA 
as a tool to maximize it is not essential here, and of course, one can use the 
different methods such as simulated annealing), we maximized it under the 
'zero-collisions' and 'no breaking-up' constraints. In this sense, this simple 
maximization process under two constraints leads to the empirical value of 7. 
Therefore, we might conclude that the flock is formed so as to maximize the 
7-value and to satisfy the two constraints. 

In fact, it is obvious from the definition that the upper bound of the 7-value 
is 1. However, from the result shown in Table. 3, the 7-value observed in our 
simulations is lower than the bound, namely, 7 ~ 0.8. We also examined to 
what extent the 7-value increases when we does not require any 'zero-collision' 
constraint during the GA dynamics and found that the 7 value increases up to 
near the bound. From these results we may conclude that the 'zero-collision' 
constraint reduces the 7-value considerably. However, we should notice that 
7-value calculated in the empirical findings by Ballerini et al. [14] also takes 
the value around 0.8 which is very close to our result. Therefore, we might 
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Generation 

Fig. 14. Evolution of 7-value in generations for the topological model. 

conclude that finding the optimal weights for the interaction by maximizing 
the 7-value under 'zero-collision' constraint is a reasonable way to realize more 
'realistic' flocking simulations even in our personal computers. 

Mathematically, inspired by the so-called Gardner's capacity in the research 
field of neural networks [20] , we might examine the following fraction of solu- 
tion space: 



1/(7) = 



S~dje(\j 


- j)5(N c )5(N h )5^(J) - 7) 


Io°°dje(\j 


-3) 



(38) 
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where €)(•••) and £(•••) stand for a step function and a delta function, respec- 
tively. Here we also defined dJ = dJid^dJ^, 0(| J\ — j) = 0(| J\ \ — ji)©(| J2I — 
j2)©(|^3| — js)- N c and AT& mean the numbers of collisions and breaking-up. 
Fixed constant variables j k) k = 1,2,3 specify the supports for these three 
variables J k) k = 1, 2, 3. The fraction v might shrink (perhaps, monotonically) 
to zero when we increases the 7 and a 'non-trivial' theoretical upper bound j c 
might be obtained as a solution of ^(7 C ) = 0. In our forthcoming article, the 
details of this argument and the results will be reported. 

11.3 Symmetry breaking in the space of interactions 

In our BOIDS modelling, we gave the interactions such as Ji, J 2 and J 3 for 
the metric model ( Ji, J 2) J3 and J4 for the topological model) to all agents as 
the same values. However, of course, one could modify the situation so as to 
give the 'agent-dependent interactions' as J% \ J^) ^° the system. Then, 
it is very interesting for us to make two-dimensional histograms for each of 
the interactions like Ji(0, 0), J<i(§, 0), ^3(0, </>)• From the histograms, we might 
obtain the useful information about the correlation (relationship) between 
anisotropy in position and anisotropy in the interaction. However, apparently 
the number of parameters to be determined by GA increases from three to 
3N and the number (system size) is critical for our computational cost of 
determination by GA within a reliable precision and realistic computational 
time. Therefore, we would like to address this problem in our future studies. 

11.4 Comparison with empirical findings 

From Fig. 8 and the same plot given by Ballerini et. al [14], 7- values for n = 1 
are almost the same and the both decrease as n increases. Moreover, these 
two curves (ours and Ballerini's) converges to 1/3 in the asymptotic limit of 
n —> 00. However, in our evaluation, the 7- value becomes lower than 1/3 in 
the range of 7 < n < 14, whereas the empirical evidence indicates 7-value 
monotonically decreases as n increases and is saturated to 1/3 beyond n = 6. 
This is an essential difference between the Ballerini's empirical findings and 
ours. 

Here, we assume that the metric interaction might cause a regular-polygon 
structure (a locally crystallized structure) in the flocking leading up to the 
'anti- anisotropy'. Hence, in this paper, we also carried out the BOIDS simu- 
lations in which the neighbours for an arbitrary agent is defined topologically 
instead of metrically. We found that the anisotropy measurement does not 
show the anti-anisotropy and the gap between our modelling and empirical 
findings was actually reduced. However, the 7-value for n — 1 is apparently 
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smaller than the empirical evidence (see Fig. 12 (left)) and it might be needed 
to introduce different types of constraints to design the artificial flocking in 
computers. 



12 Summary 



In this paper, we utilized GAs to maximize the 7- value in order to deter- 
mine the weights for interactions in the BOIDS under 'zero-collision' and 'no- 
breaking-up' constraints. We found that this procedure enables us to simu- 
late the realistic flocking phenomena even in our personal computer level. We 
showed that the resultant 7-value as a function of n-th order of the nearest 
neighbouring agents is quite similar to the empirical findings in several as- 
pects [14,21]. We carried out the simulations for both metric and topological 
definitions of neighbours in the flocking and found that the topological model 
reproduce the physical quantities of empirical evidence much better than the 
metric model does. Of course, there still exists a gap between the empiri- 
cal evidence and our results by the BOIDS modelling. Recently, Bode et. al 
[23] measured this anisotropy in artificial flockings by using stochastic, asyn- 
chronous updating scheme for each agent's movement. Their model partially 
succeeded in reproducing the behaviour of the empirical findings and showed 
us a possibility that there exists more realistic flock simulations than the con- 
ventional deterministic one. In our modelling, such probabilistic ingredients 
were not taken into account, however, their numerical evidence might stress 
that such a stochastic agents is one of essential factors to generate much more 
realistic collective behaviour of flockings. The extensive studies to reduce the 
gap should be needed to reveal the nature of non-trivial collective behaviour 
in flocking phenomena. 
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Appendix 

A Details of our genetic algorithm 

We here explain our GA procedure in this paper as follows. 

1. Initialization: Create gene configurations U. 

2. Repeat the following procedure (g — 1 : g < G). 

(a) Crossover: To make a crossover to generate a new gene configuration 
E. 

(b) Mutation: Define a gene configuration I = U + E and pick up 
| / s | + |i a | components randomly from i and mutate them. 

(c) Selection: Select the gene configurations having high fitness values 
from i. We select such \U\ gene configurations and update the pre- 
vious set U . 

We next show the above Initialization more precisely as follows. 
Initialization: 

(11) We give a random variable in the range [0.001, 0.999] to each gene Ji, J 2 
and J 3 and generate a gene configuration. 

(12) For the generated gene configuration, we evaluate the 7-value. 

(13) If the 7-value is no more than 7^ x and any 'collision' or 'breaking-up' 
is not observed, we choose the 7-value as the fitness and add the corre- 
sponding gene configuration (J 1? J 2 , J3) to the set U. 

(14) Repeat the above (il),(i2) and (i3) until the number of configurations 
reaches \U\. 

Where we defined the 7^] x to confirm that one can get optimal weights even 
if he (or her) starts the GA from wrong initial gene configurations having 
relatively low fitness values. We should keep in mind that the evaluation of 
7-value is carried out for 8-independent runs of simulations and we set the 
7-value which is averaged over 8-independent runs to the fitness function if 
and only if there is no 'collision' or 'breaking-up' during each trial. 

We next explain the details of the Crossover. 
Crossover: 

(cl) Pick up arbitrary two gene configurations from the set U and define these 

configurations as a and b. 
(c2) We swap arbitrary one gene in the a for arbitrary one gene in the b. 
(c3) We evaluate the 7-value for the modified a, b and add them to the set E 
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if and only if there is no 'collision' or 'breaking-up'. 
(c4) Repeat the above procedures (cl),(c2) and (c3) until we have \E\ gene 
configurations. 

As the searching (solution) space is constructed by only three variables Ji, J 2 
and J 3 having real numbers, it is naturally assumed that the effect of the 
crossover is relatively weak on the optimization by GAs. Thus, we overcome 
this weakness by introducing the following two kinds of Mutations. 

Mutation 

(ml) From the set / = U + E, we pick up a gene configuration and define it 
as c s . 

(m2) We update arbitrary one gene, say, J i) in the configuration c s randomly as 
Ji —> Ji+S, where S stands for a random number in the range [0.001, 0.01]. 

(m3) For the modified c 5 , we evaluate the 7- value and add the c s to the set / 
if and only if there is no 'collision' or 'breaking-up'. 

(m4) Repeat the above (ml),(m2) and (m3) until we obtain |/ s | gene configu- 
rations. 

(m5) We pick up an arbitrary gene configuration and replace a single gene Ji 
with a random number in the range [0.001,0.999]. We make |/ a | gene 
configurations by making use of the same procedure and add the set 7 a 
to the set /. 

It should be noted that the above procedures (ml),(m2) and (m3) achieves 
'local search' whereas the procedure (m5) acts as 'global search'. Thus, the 
above Mutation realizes the effective searching by using a mixture of local 
and global searches. 

Finally, we explain the details of the Selection as follows. 
Selection 

(si) Pick up a gene configuration having the highest 7- value from the set I as 
an elite and the evaluate the 7-value if and only if there is no 'collision' 
or 'breaking-up'. Then, we add the gene configuration to the set U in the 
next generation. 

(s2) For all genes in the set /, we calculate the 7-values as the corresponding 
fitness values. If g > G/2, we make a linear transform on the 7-value 
for all gene configurations in the set I and choose the transformed set of 
7-values as fitness functions. 

(s3) Make a roulette selection of the gene configuration based on the fitness 
functions and evaluate the 7-value if and only if there is no 'collision' 
or 'breaking-up'. Then, we add the gene configuration to the set U and 
repeat the above procedure until we have \U\ configurations. 
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Number of generation (G) 


25 


Size of gene configurations (\U\) 


1U 


T"n 1 "1" 1 CI 1 0/-TY1flY 

niiLidi y iiidA 




Crossover rate (\E\) 


i 


Mutation rate 1 (|/ s |) 


2 


Mutation rate 2 ( 7 a ) 


1 



Table A.l 

Parameter setting for our GA. 

As we are restricted ourselves to the case in which the number of trials, the 
time of observation are limited, there exists a possibility that we get un- 
expected high 7-value and the gene configuration having such high 7-value 
accidentally and to make matter worse, such gene configuration sometimes 
might survive until the final generation. To overcome this difficulty, we re- 
peat the measurement of the 7-value for the selected gene configurations until 
'zero-collision' condition is strictly satisfied and we delete the selected gene 
configurations if they lead to 'collision' or c breaking-up' even if they show the 
high 7-value. Thus, we systematically solve the optimization problem under 
two essential constraints. 

A.l The size of GA in computer simulations 

We next explain the size of GA in our computer simulations. We summarize 
them in Table A.l. We carry out the GA having the above setting-up of the 
parameters. Then, for each generation of the GA, we evaluate the 7-value, the 
weight of the interactions and the distribution of the gene configurations for 
three independent runs. 



B Border bias effect and procedure to avoid it 

In computer simulation for finite size systems (N <C 00), we should keep in 
mind that the results are sometimes influenced by the so-called border bias ef- 
fect [22]. The border bias effect comes from the asymmetric shape of flockings. 
For instance, the asymmetric shape we mentioned here is typically an elliptical 
shape which is obtained as a deviation from a given symmetric sphere. The 
major axis of the elliptic corresponds to the direction of moving of the flock. 
For a given flocking and a given agent, the n-th nearest neighbouring agent is 
more likely to exist in the direction of moving of the flocking rather than the 
direction perpendicular to the flock's movement. As the result, the 7-value 
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decreases to below 1/3 as n increases, and then, 'anti-anisotropy' emerges. 
This phenomenon is nothing but border bias effect we mentioned here. It is 
very important for us to avoid the border bias effect to evaluate the 7-value 
precisely. If there is no correlation between the direction of flock's moving and 
the shape of the flocking, one may cancel the effect by taking the average of 
the 7-value over several (usually, a lot of ) trials with different initial condi- 
tions. However, if not, it is very difficult for us to cancel the effect by this 
simple procedure. In most cases of simulations, the shape of flockings is an 
elliptic in which the major axis is always in the direction of flock's moving and 
there exists an apparent correlation between the shape and the direction of 
moving. In general, the effect is more serious for the flocking simulation with 
small number of agents than that with huge number of agents. Ballerini and 
Cavagna et. al. cancelled the effect by excluding the agents on the border, 
however, the usefulness of their procedure is limited to the case in which the 
number of mates is larger than 400 [14,22]. 

From this fact in mind, we provide two distinct procedures to cancel the two 
types of border bias, namely, multiple trial method for border bias caused by 
'cubic shape' (slab) of the flocking and sphere extraction method for border 
bias induced by 'spherical shape' of the flock, we show that these procedures 
make our flock simulations free from the border bias effects. 

B.l Multiple trial method 

We first examine the extreme situation which shows 'fake '-anisotropy by bor- 
der bias [22]. 

Let us think about the artificial flocking having a cubic shape as shown in the 
left panel of Fig.B.l. The ratio of three slides of the cube is given as 7 : 3 : 1. 
Then, we evaluate the 7-value for the flock moving to the direction of the 
shortest side. We also purposely make the flocking so as to show the isotropy 
for the angular distribution for individual by hand. Hence, we should naturally 
observe the 'isotropy' through the 7-value, namely, 7 = 1/3 if one correctly 
simulates the artificial flocking without any border bias effect. 

The result is shown in the right panel of Fig. B.l by 'circles'. From this panel, 
we clearly find that the anisotropy emerges and our simulations are apparently 
affected by the border bias. Even if we change the direction of flock's motion 
from the shortest side to the longest slide, one cannot obtain the 'isotropy' 
but one observes the 'anti- anisotropy' as shown in the right panel of Fig. B.l 
by 'boxes'. 

In order to overcome this type of border bias effects, we utilize the so-called 
multiple trial method. Namely, we evaluate the measurements such as the 7- 
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Direction given by the random number 



2 4 6 8 10 12 14 
order of the neighbour (n) 

Fig. B.l. The extreme shape ('cubic shape') of artificial flocking which shows 
'fake'-anisotropy (left). The right panel shows the 7- value as a function of order 
n. Circles and boxes correspond to the 7-values observed in the flocking moving 
to the directions of the shortest and the longest sides, respectively. We purposely 
put 1200-individuals randomly into this cubic so as to show the isotropy 7 = 1/3 
by hand, however, the anisotropy/'anti-anisotropy' emerges. The 'diamonds' in this 
panel show the result by the multiple trial method. We clearly find that the isotropy 
is observed as we expected. 

value for a lot of independent trials, in each of which we give the initial di- 
rection of flock's motion randomly without any correlations with any specific 
direction such as the shortest and the longest directions of cubic. Then, the 
measurement is obtained as an average over the trials. We show the resulting 
7-value for 1000-trials in the right panel of Fig. B.l by 'diamonds'. We clearly 
find that the 'isotropy' is observed as we expected. 

B. 2 Sphere extraction method 

In the previous subsection B.l, we were confirmed that the multiple trial 
method efficiently reduces the border bias effect. However, we should notice 
that the method is effective only for the case in which the direction of flock's 
movement (specified by the velocity vector V) and the shape of the flocking 
has no correlation. 

To understand it, let us define a unit vector pointing to the direction of the 
longest side of the cubic in the previous subsection B.l by e^. Then, the 
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Fig. B.2. The 7-value as a function of order n. A typical snapshot of the aggregation 
is shown in the inset. 



following condition should be satisfied 



e L - V = (B.l) 

to use the multiple trial method effectively, where (• • • ) stands for the time 
average over the sampling from BOIDS dynamics. It is obviously confirmed 
when we consider the extreme case in which randomly given vectors Vk (k = 
1, • • • , M: M is the number of trials) as the direction of flock's motion are 
always correlated with in terms of • Vk — 1 {k — 1, • • • , M). For this 
case, we clearly notice that the multiple trial method gives exactly the same 
result as the 'circles' shown in the right panel of Fig. B.l. 

Therefore, another way to reduce the border bias effect is needed to evaluate 
the anisotropy measurement precisely in our BOIDS simulations. 

In order to reduce the border bias effect in such cases, we take only agents who 
are in the sphere being inscribed with the shape of the flocking to calculate 
the 7-value. To check the usefulness of our procedure, we set-up 100 agents 
moving with the same velocity in the same direction and give them their 
initial positions uniformly in the elliptic in which the major axis is in the 
direction of flock's moving (see the inset of Fig. B.2). For these agents, we 
evaluate the 7-value up to n = 25 and compare the result with the 7-value 
calculated without the above procedure. The results are shown in Fig. B.2. We 
carried out 1000-independent runs and calculated the average of the 7-value. 



39 



As we mentioned, the distribution of the position of agents is uniform and 
the correct (true) 7- value should be 1/3. From this figure, we find that our 
procedure (circle Sphere) works well in comparison with the result without the 
procedure (diamond Normal). 
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